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Abstract 

The synchrosqueezing method aims at decomposing ID functions as superpositions 
of a small number of "Intrinsic Modes", supposed to be well separated both in time 
and frequency. Based on the unidimensional wavelet transform and its reconstruction 
properties, the synchrosqueezing transform provides a powerful representation of multi- 
component signals in the time-frequency plane, together with a reconstruction of each 
mode. 

In this paper, a bidimensional version of the synchrosqueezing transform is defined, 
by considering a well-adapted extension of the concept of analytic signal to images: the 
monogenic signal. The natural bidimensional counterpart of the notion of Intrinsic Mode 
is then the concept of "Intrinsic Monogenic Mode" that we define. Thereafter, we inves- 
tigate the properties of its associated Monogenic Wavelet Decomposition. This leads to a 
. natural bivariate extension of the Synchrosqueezed Wavelet Transform, for decomposing 

and processing multicomponent images. Numerical tests validate the effectiveness of the 
method for different examples. 



| Keywords: monogenic signal, wavelet transform, directional time-frequency image analysis, 

synchrosqueezing . 
MSC 2010: 65T60, 92C55, 94A08. 
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1 Introduction 

In the past few years, there has been an increasing interest in representing images using 
spatially varying sinusoidal waves. As the understanding of the theory advanced, amplitude- 
and frequency-modulation (AM-FM) decompositions have been applied in a large range of 
problems: for example motion estimation using a flow-optical method based on an assumption 
of phase-invariance [H [3] , reconstruction of breast cancer images [U [5] , texture analysis [6j 
E], or ultrasound image segmentation methods [8] based on the concept of monogenic signal 
and quadrature filters [SJ. A survey of applications of AM-FM decompositions in medical 
imaging can also be found in [lOj . 

In each case, the main challenge is to decompose any input image s(x) into a sum of 
bidimensional AM-FM harmonics of the form 

L L 

(1) s(x 1 ,x 2 ) = y^^t) = } y Aj(x 1 ,x 2 ) cos(y»f (g?i, x 2 )) , 

i=i i=i 
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where An > denotes a slowly-varying amplitude function, tpi denotes the phase, and t = 
1, • • • ,L indexes the different AM-FM harmonics. To each phase function, one can associate 
an instantaneous frequency vector field defined as uj£ = V<pg. Finding the components S£ from 
the bidimensional signal s is called the decomposition problem. 

When the decomposition is given, another subsequent problem is to determine for each 
component its amplitude, phase and frequency functions involved in ([1]). The amplitudes 
can be related to the energy contained at each point of the image, whereas significant texture 
variations are captured in the frequency content. For example, in the case of a single compo- 
nent, the instantaneous frequency vectors are orthogonal to the isovalue intensity lines of an 
image, while their magnitudes provide a measure of local frequency content. The problem of 
estimating the amplitudes, phases and instantaneous frequencies of a signal of the form (pQ) 
(bidimensional or not) is called the demodulation problem. 

In the one dimensional case, these two problems have been widely investigated both from 
the theoretical and practical point of view. To address the demodulation problem in the one 
dimensional setting, a first and necessary step is to define, in a proper way, the concept of 
amplitude, phase and instantaneous frequency of a given signal. To this end, one can use the 
well-known Hilbert transform, defined for any / 6 L 2 (R) as 

For a.e. t G M, Hf(t) = lim ( - [ ds 

e^O y-K J\t- S \ >£ t - S 

The analytic signal associated to / is then the complex- valued function F = {1 + \H)f . 
Thereafter, one defines amplitude A and phase (p of / as, respectively, the modulus and 
argument of the analytic signal F associated to / (the uniqueness of phase being ensured 
under smoothness assumptions on / and under some initial conditions of the form <p(0) = ipo). 
One then obtains 

For a.e. t 6 R, F(t) = A{t)e lLp{t) . 

The instantaneous frequency of / is thus the derivative of the phase ip' [llj . The practical 
estimation of amplitude, phase and instantaneous frequency was the subject of numerous 
studies (see [12] for a survey). The "naive" method, which consists in using the Hilbert trans- 
form to estimate A and tp, is known to be numerically unstable [13]. Alternative methods use 
the wavelet transform associated to reallocation techniques: in the one dimensional case, the 
time-frequency localization of wavelets allows to compute robust estimations of instantaneous 
frequency lines in the time-frequency representation: several approaches have been developed 
in the 90th, one should cite the method of "wavelet ridges" (see the two pioneer works [14J 
and [H>]), the squeezing method introduced in |16J or the reassignment method [17j . 

On the other side, the decomposition problem consists in developing decompositions of 
any signal into a sum of "well-behaved" AM-FM components, and was widely studied. A 
recent attempt in this context is Empirical Mode Decomposition (EMD), initially introduced 
by Huang et al. in [18] and later popularized by Flandrin and his co-authors (see e.g. [19] 
or [20]). 

Basically, EMD is the output of an iterative algorithm which provides an adaptive decom- 
position of any signal into several AM-FM components. It is now used in a wide range of 
applications including meteorology, structural stability analysis, and medical studies [T2J. In 
spite of its simplicity and its efficiency, this method is hard to analyze mathematically, since 
it is defined in an empirical way. To tackle this problem, in [21], the authors propose and 
analyze an alternative method called "synchrosqueezing" derived from reassignment methods: 
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they introduce a class of functions which can be viewed as the superposition of a reasonable 
number of modes, i.e. which read Yle Ae(t)e ilpe ^ , where the amplitude Ag(t) of each mode 
is slowly varying with respect to its instantaneous frequency wg(i) = (p'At). For each func- 
tion belonging to this class, the synchrosqueezing provides an estimation of this (unknown) 
decomposition. Synchrosqueezing presents many similarities with Empirical Mode Decompo- 
sition as shown in [21] (see also [22J for a comparison of these two methods or [23] for more 
about potential applications). In addition, since synchrosqueezing is based on the wavelet 
transform, it can be used to address simultaneously the decomposition and demodulation 
problem as done in |21j . 

In this paper, we will focus on the two dimensional case and extend the synchrosqueezing 
approach to images. Our first step consists in defining a convenient extension of the concept of 
analytic signal, in order to properly define the notion of amplitude and phase of an image. Two 
main generalizations of the analytic signal to the bidimensional setting have been considered in 
the literature: the hypercomplex and the monogenic signal defined respectively in [21] and [9]. 
Here, we focus on the monogenic setting, since its characteristics lead to easier interpretation. 
The monogenic signal associated to an image is defined using the Riesz transform. If we are 
given a real valued function / 6 L 2 (R 2 ), one defines its Riesz transform IZf as the following 
vector-valued function 



where for any i = 1, 2 and a.e. i£l 2 one has 




'£<■/(■'■) lim . I - I rrl/fe) dy 

\x~ y\>e 1*^ y\ j 



Then, the monogenic signal associated to / is 

Mf = 



f 

Kf 



The characteristics of the monogenic signal are usually defined using the quaternionic formal- 
ism (see Appendix [A] for some background about quaternionic calculus and Appendix [B] for 
the main properties of the Riesz transform): the monogenic signal Aif associated to / can 
be read 

Mf = f + n 1 fi + n 2 fj, 

where (1, i, j, k) denotes the canonical basis of the algebra H of the quaternions. In [9j and [25J, 
it is proved that the monogenic signal can be written as 

Mf = A e^™ 8 = A (cos (p + ng sin ip) with ng = cos 9 i + sin 9 j , 

where A is a positive real valued function and <p and 9 are two real valued functions defined 
in a unique way (under suitable assumptions). The function A is called the amplitude of /, (p 
its phase and ng its local orientation. The function uj = Vip will be called the instantaneous 
frequency of /. To illustrate these definitions, we consider the case where f(x) = Aq cos(k ■ x) 
with k = (ki, ^2) an d k\ > 0. Set 9q = Arctan(A;2/A:i). One then has 

\ a fsin(k 1 x 1 + k 2 x 2 ) cos 6> \ . k , . 
\ v sm(A;iXi + ft 2 x 2 ) smfc'oy \k\ 
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and 




— A e (k-x)(cos0 i+sin6> j) 



Then for all x G M 2 



one can set 



A(x) = Aq, <p(x) = k ■ x, 9{x) = 9q = Arctan(/c2/^i) • 



The aim of the present paper is to define the synchrosqueezing of bidimensional images, 
in the monogenic setting, in order to address both the decomposition and the demodula- 
tion problems. To this end, we will consider monogenic signals which are superposition of 
monogenic waves of the form 



where A and 9 are slowly varying functions with respect to ip. 

We will consider successively the two different problems of deconvolution and decompo- 
sition. As in [21j, our method is based on the continuous wavelet transform. Therefore, in 
Section [2J we recall a recent extension of the continuous bidimensional wavelet transform 
to the monogenic setting. We then investigate in Section [3J the problem of estimating the 
amplitude and instantaneous frequency of a monogenic spectral line of the form ([2]). Subse- 
quently, we define in Section H] the bidimensional synchrosqueezing and its application to the 
decomposition of multicomponent images. Finally, Section [5] gives some insights on the prac- 
tical implementation of bidimensional synchrosqueezing, and provides numerical experiments. 
Proofs are postponed to Section [6] to make reading easier. 

2 2D-wavelet transform and monogenic wavelet transform 

In the one dimensional setting, the synchrosqueezing method is based on the one-dimensional 
wavelet analysis: indeed, the wavelet transform is used both to decompose an analytic signal 
into several components, and also to estimate the amplitude and the instantaneous frequency 
of each component. The main advantage of using the wavelet analysis is to avoid the use of 
the discrete Hilbert transform, since there is a link between the wavelet coefficients of the 
analytic signal associated to a real function /, and its analytic wavelet coefficients (i.e. the 
wavelet coefficients of / among an analytic wavelet). 

In the bidimensional case, it is then quite natural to wonder if this approach can be 
extended using bidimensional wavelet analysis, replacing the analytic signal by the monogenic 
signal associated to our data. Since the monogenic signal can be viewed as a 3D vector field, 
it can be componentwise analyzed by the usual wavelet transform: Section 12.11 recalls first 
some basics about the bidimensional continuous wavelet transform. Thereafter in Section \2.2\ 
we recall the main characteristics of the monogenic wavelet analysis, which is an extension of 
the usual wavelet analysis introduced in [26] and |27| . We would like to point out that, for 
any bidimensional real image /, the wavelet coefficients of the monogenic signal Mf can be 
related to the monogenic wavelet coefficients of / (see Proposition 12 . 1 1 for a precise statement). 
Then, using monogenic wavelet analysis, one can recover the characteristics of the monogenic 
signal associated to our data, without estimating it. 

In what follows, for any (a, 6, a) 6 x R x (0, 2ir), we denote by D a the dilation operator, 
T;, the translation operator and R a the rotation operator defined on L 2 (R 2 ) by: 



(2) 




D a f(x) 



a- 1 fix/a), T b f(x) = f(x - b), R a f(x) = f{r~ x x) 
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where r a is the usual 2x2 rotation matrix of angle a: 

/cos a —sin a 

(3) r a = 

ysm a cos a 

2.1 The usual bidimensional wavelet transform 

We briefly recall classical definitions about the bidimensional continuous wavelet transform. 
A (real or complex) function ip G L 2 (M. 2 ) is called an admissible wavelet if it satisfies 0: 

(4) C^ = (2nf r ifigfu < +oc 

As usual (see [28]), the wavelet family {^ a ,^&}( a ,a,6)eR* x(o,2tt)x]R2 is defined by ip a ^ b = 
T b R a D a i\). One then defines the wavelet coefficients of any / G L 2 (M 2 ) as follows: 



c f (a,a,b) = / f(x) tp a ,a,b(x) dx . 
Jk 2 

If the wavelet is assumed to be isotropic, the wavelet coefficients do not depend on a. In this 
case, we will denote c/(a, b) = c/(a, 0, b) and Tp at b = ip a ,o,b- 

Any square integrable function / can be recovered from its wavelet coefficients using the 
so-called reconstruction formula: 

(5) f(x) = ^- f f f c f (a,a,b) ip a>a ,b(x) ^da db , 

U V ^6eR 2 iaG(0,27r) JaG(0,+oo) a 

where this equality stands in L 2 (M 2 ). When the wavelet is assumed to be isotropic, we get a 
simpler expression: 

f( x ) = TT [ [ c f (a,b) ip a ,b(x) ^dfe . 

W JbeK 2 JaG(0,+oo) a 

Other reconstruction formulas are available (see e.g. [28]), among which the pointwise recon- 
struction formula, obtained by summing over scales only. In the isotropic case, it reads: 



2?r f +co . . da . ~ f 



(6) f(x) = 7 r- c f (a,x) -5- with = / d£ . 

Jo a Jr 2 I? I 

As in the one dimensional case [21], this formula will play a special role in the bidimensional 
synchrosqueezing. 

For any real function / G L 2 (M 2 ), one can also define the wavelet transform of its monogenic 
signal F = Mf = f + K x f i + K 2 f j , as: 

■ ( c < 

c F = c f + c nif 1 + cn 2f j = I c nif 

W2// 



lr rhe Fourier Transform of if> being defined by: V>(£) = ^ J R2 tp{x)e '^' x dx 
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2.2 The monogenic wavelet transform 

We now present the continuous monogenic wavelet transform as denned in We consider 

(*\ 

a real admissible wavelet ip and define ip( M > = A4i[) = I TZiip I the associated monogenic 

wavelet. Observe that %fj^ M ^ is also an admissible wavelet (taking values in R 3 ), since each of 
its components satisfies relation (j4]). We define below the monogenic wavelet coefficients of a 
real function. 

Definition 2.1 Let f £ L 2 (R 2 ). The monogenic wavelet coefficients Cj^^a, 6) of f are 
defined by: 

(7) cf 4) (a,a,b)= [ f(x)^ a M J b (x)dx, 

where for any (a, a, b) <G (0, +00) x (0, 2vr) x M 2 , = T b R a D a (Mip). 

Remark 2.1 For any (a, a, 6), c^\a,a,b) is a Clifford vector, and thus can be written as 
Cj M ^(a, a, b) = Cf(a, a, b) + c'p (a, a,b) i + c'p (a, a, b) j where we denote for i = 1, 2 

cy(a,a,b) = / /(x) (T b R a D a Kiip)(x) dx , 

('see Appendix\M for more details). 

The next proposition states that there exists an explicit relationship between the mono- 
genic coefficients of a real image / and the usual wavelet coefficients of its monogenic signal 
Mf. 

Proposition 2.1 Let f £ L 2 (M 2 ) and denote by F = Aif the monogenic signal associated 
to f. Then for any (a, a, b) G (0, +00) x (0, 2ir) x M 2 



(8) cf(ci, a, b) 

Proof. We will prove that 

C/ M) (°' a ' 6 ) = (J Ci? ( a ' a ' 6 ) ' 

which is equivalent to (|8j). By definition of the monogenic wavelet transform one has: 

c ( f M) (a,a,b) = [ f{x) T b R a D a (MiP)(x) dx . 

The translation-, scale-invariance and steerability properties of the Riesz transform recalled 
in Propositions IB.31 and IB~4l of Appendix [Bl imply that: 

T h R a D a {n^) = r- l n{^ b ) . 
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Hence 
(9) 



f(x) T b R a D a (K*p)(x) dx = r~ l / f(x) n{i> a ^ b ){x) dx . 



Since by Proposition IB. 51 the Riesz transform is a componentwise antisymmetric operator on 
L 2 (R 2 ), one deduces that: 



(Kf)(x) 1p a ,a,b{x) dx . 



Equations ([9]) and ([10]) then directly imply relation 



□ 



We now consider the case an admissible isotropic real-valued wavelet ip. Examples include 
the Mexican hat (Laplacian of Gaussian) or the Morse wavelets (see |29j and Section [5]) . In 
such a case, the wavelets tp a ,a,b as well as the wavelet coefficients cf(cl, a, b) of the monogenic 
signal do not depend on the orientation a. From now, we then denote tp a ,a,b and cp(a,a,b) 
respectively by ^> 0) & and cp(a,b). Thus equation ([8]) can be simplified in the following way: 



c F (a,b) 





-r, 



a, a, b) 



c f (a,b) 





-TV 



/ c f (a,b) 

Cj\a, q, b) 
\cf\a,a,b) 



- cos a ci (a, a, b) + sin a (a, a, 6) 
y— sin a Cy. (a, a, 6) — cos a ci 2 (a, a, 6) / 

In practical situations, it will be easier to consider the case in which a = 0. We then get that 



(12) 



c F (a,b) 







c ( , M) (a,0,6). 



To sum up, in the isotropic case, the monogenic wavelet coefficients of the 2D real image / can 
be related in a very simple way to the wavelet coefficients of the monogenic signal using f)12|) . 
It will be useful in the sequel since when analyzing images, one deals with a real-valued signal 
/ and aims at recovering the wavelet coefficients cp(a,6) of its monogenic signal F without 
computing F. 



3 Wavelet analysis of bidimensional spectral lines 

In this section, we tackle the problem of estimating the amplitude and instantaneous fre- 
quency of a bidimensional mode of the form A(x)e !fi ^ n ^ x \ In the one dimensional setting, the 
"wavelet-ridge" method has proved to be efficient to address this problem |14tll5j. In |26tl29| 
it has been extended to the bidimensional setting using monogenic wavelet analysis. 

Here we follow an alternative approach consisting in extending synchrosqueezing to the 
bidimensional context. We first define precisely the bidimensional modes that we are consid- 
ering, and that we will call "Intrinsic Monogenic Mode Function" (IMMF). Then we deduce 
the notion of bidimensional frequency vector and state an approximation result using the 
monogenic wavelet transform. 

Let us first define our concept of bidimensional modes: 
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Definition 3.1 Let e > 0. An Intrinsic Monogenic Mode Function (IMMF) with accuracy e 
and smoothness a > is a function F G B^ T l (IR 2 , M) of the form 

(13) F(x) = A(x)e^ x)nB ^ with n e(x) = cos(0(x))i + sin(#(x))j , 

where A>0andA,6£ C 1 (M 2 ) n L°° (R 2 ) , <p G C 2 (IR 2 ). The function A is called the amplitude 
of F, whereas (p and n$ are called respectively the scalar phase and the local orientation of F. 
The function ip is assumed to have bounded derivatives of order 1 and 2: 

(14) Vi G {1,2}, < inf \d Xi <p(x)\ < sup \d Xi <p(x)\ < oo . 



< oo . 



(15) M = max sup d 2 . x . ip(x 

Further we assume that A, 9, X7 2 tp are slowly varying functions with respect to V(p, namely 
that for i = 1 , 2 

(16) \fx G R 2 ,max(\d Xt A(x)\,\d Xt 9(x)\) < e\d Xi <p(x)\ . 
and 



(17) Vx G M 2 , . max n d*M x ) < ■ 



ii,«2=li2 



)2 



Remark 3.1 Observe that, in the sequel, we shall consider non-square integrable functions. 
In this case, the two reconstruction formulas |2P and may not hold (see Appendix B 
of }30^). We then need an additive assumption to ensure that the integral f + °° cf(o, 6)^ 
exists to give some sense to the reconstruction formula Moreover, we will estimate the 
behavior of 

/ \cF{a,b)\ — 

for any ao > 0. 

The additional hypothesis F G B^ l (IR 2 , H) for a > is equivalent to the following condition 
on the wavelet coefficients (see for e.g. ]31^): 

f°° da 

(18) IIFIId-o- = / sup I cf (a, b)\—z — < +oo . 

°°> l Jo bm 2 a a 

In such case, we easily deduce that for any ao > 0: 

f°° da f°° 1 da 

(19) sup/ \c F (a,b)\^ = sup / \c F (a, b)\ — ^— < a a \\F\\^-a 
bm 2 Ja a bm 2 Ja a* a 2 ° a oo,i 

Moreover, by hypothesis, F G C 1 (M 2 ,EI) and VF is bounded. Using the expression of the 
wavelet coefficients then implies that: 

\c F (a, b)\ < a 2 ||VF|| LO c [|a^|| £ i , V6 G IR 2 

which ensures that for any oq > 

f a ° da 

(20) sup / \c F (a, b)\— < oo . 

beR 2 Jo a 

and finally leads to the existence of the integral f + °° cp(a, b)^. 
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Remark 3.2 In ]29^ authors considered spectral lines of the form f{x) = A(x) cos((p(x)) with 
A,Vip slowly varying with respect to (p. The authors then approximate the monogenic signal 
F associated to f by an IMMF of the form 1113]) with constant amplitude, linear phase and 
constant local orientation. 



We now focus on the estimation of the amplitude and instantaneous frequency of an IMMF. 
To this end, we need some assumptions on the wavelet that we list below: 
Assumptions (W) 

(i) the wavelet ip is isotropic and real-valued. 

(ii) V e W 1 ' 1 (R Z ). 

(hi) The hrst three moments of tp and are hnite, namely 

sup I a < oo, sup l' a < oo , 

aG{l,2,3} ae{l,2,3} 

where for any a > 

(21) I a = [ |x| Q |V>(x)|dx , I' a = f \x\ a \Vip(x)\dx . 

JR2 JR2 



We first deal with the simple case of an IMMF of constant amplitude and phase: 

Lemma 3.1 Let F be an IMMF of the form F(x) = ^ e (fc ' x)ns with k = (ki,k 2 ) G R 2 , 
Aq G R+,# G R, and let ip be an isotropic real wavelet belonging to iy 1 ' 1 (R 2 ). 
The wavelet transform of F is given by: 



(22) c F (a,b) = A a^(ak)(cos(k-b) + sm(k-b)(cos9 i + sin#j)) = a^(ak) (A Q e {k - b)ne 
and one has for i = 1,2: 

(23) d h c F (a, b) = hn e (af (a*;)) (A e^ n ^ . 

The vectors k and ng can then be recovered using the two following relations: 

kin e = d bl c F (a,b) x (c^a,^) -1 , 
k2ng = d b2 CF(a,b) x (cj?(a,6)) _1 . 



Proof. By definition, since the wavelet is isotropic real, one has: 

c F (a, b) = a -1 / F(x)ip ( ] dx = a f F(au + b)ip(u)du 



A simple calculation leads to: 

(J R2 cos(k ■ (au + b))ij}(u)du \ / cos(k.b) 

cos 9 f R2 sin(fc • (au + b))ij}(u)du = AQa^(ak) sin(k.b) cos(9) = aip(ak)F(b) . 
sin9 L a sin(k ■ (au + b))ip(u)du J \sm(k.b) sin(0) / 
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where we have used ip(-ak) = ip{ak) since the wavelet is isotropic. 
Now, since tp G W 1 ' 1 ^ 2 ) and F G L°°(IR 2 ), we get: 

d h c F {a,b) = a' 1 [ F(x)d bi L (—)] dx = -a~ 2 [ F(x)(d Xi ip) (—) dx . 

Jr-2 I V a )\ J R2 \ a J 

A similar calculation, using d Xi ip(£) = i£i*0(£),V£ = (£1,62)) and tp isotropic, yields: 

(— sin(k.b) \ 
cos(fc.6) cos(6>) = akitp(ak) ngF(b) , 
cos(fc.fr) sin(0) / 

where we have used: 

ngF(b) = ng(cos(k ■ b) + ng sin(k ■ b)) = — sin(/c ■ b) + rig cos(k ■ b). 

This leads to □ 
Remark that, more generally, if F(x) = J 4 e( fc ' a;+a ) ne , (a £ 1): 

(24) c F (a,b) = A a^{ak)e^ b+a ^ = a$(ak) F(b) 

(25) 3^0^(0,6) = akiip(ak) ngF(b) 

The formulas (|22p and (|23[) can be extended in the general case of an IMMF F of the form f)13[) . 

Proposition 3.1 iei -0 a wavelet satisfying assumptions (W) and F an IMMF with accuracy 
e and smoothness a > of the form Iil3\) . For all (a, b) G x M 2 , one has 



(i) 



c F (a, 5) = a$(aVip(b)) (A(b)e^ b)n(b) ) +ea 2 R 1 {a,b) , 



with 
(26) 



|#iM)l < Ii(V^A(6) + l)|V^(6)|+aJ 2 A(6)(|V^(6)|+V2M)/2 
+aI 2 M/2 + a 2 / 3 ^(6)M/6 . 



^c F (a,6) = d b Mb)n e(b) (a#(aVp(&))) (yl(6)e^^w) +eaii 2 (a,6) , 



with 



\R 2 (a,b)\ < A(b){a\V^(b)\r 2 /2 + a 2 MrjQ) 
{ ' + (>/2A(6) + 1) (|V^(6)|I( + aM/ 2 /2) . 

Remark 3.3 As a consequence of this proposition, since in the isotropic case c F {a, b) and 
Cj M ^ (a, a, b) are related by All]) , the above proposition specifies the results stated in \20$ 
and ]29$ . Notably, it gives a bound of the approximation error. 
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Proof. See Section EU □ 
We now define the two following Clifford vectors: 

(28) Ai(o,6) = d bl c F (a,b) x (cfM)) -1 , 

A 2 (a,6) = d b2 CF(a,b) x (cF(a,b)) L . 

These two vectors will provide an approximation of the vectors c\.(p(b)nQn,\ for an IMMF F 
of the form (|13p . More precisely, we can state the following result: 



Theorem 3.2 Let F be an IMMF with accuracy e > and smoothness a > of the 
form [T3\) . Assume that we are given tp a wavelet satisfying assumptions (W) and consider 
(a, b) G W* + x M 2 such that 

\cF(a, b)\ >e u for some v G (0, 1/2) . 

Then for i = 1, 2: 

(29) \Ai(a,b) - d h <p(b)n e{b) \ < e 1 ^ a (\R 2 (a,b)\+ a|fli(o,6)| |Vp(6)|) 



where R\, R2 have been defined respectively in [26]) and (27\ ). In addition for any oq > 0, 
there exists some Eq > such that, for all (a,b) G (0, ao) x M 2 and any < e < e^: 

(30) |Ai(a,6)-9 6 ^(6H (6) | <e v 

Proof. See Section E2J □ 



Since cos(tp(x)) = cos(— <p(x)), there is an ambiguity in the definition of the instantaneous 
frequency since ip can be replaced by — (p. To solve this problem, practitioners usually 
assume that the first component of the instantaneous frequency is positive (see [10] and 
references therein for more details). Under this additional assumption, our method yields an 
estimate of the amplitude and of the instantaneous frequency of an IMMF: 

Proposition 3.3 The notations and assumptions are those of Theorem \3.2\ In addition, we 
assume that: 

Vx G M 2 , d xi ip(x) > . 
Then for i = 1,2, and (a, b) G (0, ao) x M: 

(31) |^^(b) - |A x (a,6)|| < e" . 

and 



(32) \d b Mb) ~ |A 2 (o,6)| sgn(Re(9 fel c F (a, b) d b2 c F (a,b)))\ < e v . 

Proof. See Section E3J □ 
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4 Decomposition of a mult i— component image into spectral 
lines 



In this section we will consider a more general context, where our model is now a superpo- 
sition of several IMMFs (introduced in Section [3]), assumed to be well separated, as defined 
in the following definition 14.11 We will show in Theorem 14.11 how such functions can be de- 
composed into monogenic modes, using synchrosqueezing. Thereafter, for each component, 
the associated amplitude and instantaneous frequency will be estimated using the results of 
Section [3j 

Definition 4.1 A function F defined from R 2 to H is said to be a superposition of well 
separated Intrinsic Monogenic Mode Components with smoothness a > up to accuracy 
e > and with separation d > 0, if there exists a finite integer L such that 

L 

(33) F(x) = ^F e (x), 

where all the Fi are IMMFs with accuracy e and smoothness o~i > a > of the form i!3\): 
Fg(x) = Ag(x) Q Pe ^ n9 t( x ) ; and moreover satisfy for any x, £ > £' and i G {1, 2}; 

(34) \d Xi <pt(x)\ > \d x .tp£>(x)\ , 
and 

(35) \d Xi ipi(x)n edx) - d Xi ip e >(x)n 9i/ ( x) \ > d [\d x .<pt(x)\ + \d Xi ipt> (x)\] . 
We denote by A £ a d the class of all functions F satisfying these conditions. 

We now define the monogenic synchrosqueezing transform (MS ST) of a function F belonging 
to A Eja> d- 

Definition 4.2 Assume that we are given ip a wavelet satisfying assumptions (W) such that 

supp(V0 C {£ G R 2 ; 1 - A < |f | < 1 + A} , 

with 

(36) A < d 



2{\ + d) 



Consider h G C£°(R,R) such that f R h(x)dx = 1. Let 5 > 0, v G (0,1), and e > 0. For 
F G A e>a ,d, the monogenic synchrosqueezing transform of F is defined for any (6, k, n) G 
l 2 xl 2 x S 1 fS 1 being the unit sphere of R 2 ) by 
(37) 

6u f . 1 /fci-Re(n Ai(a,6))\ fk 2 -Re(n A 2 (a,6))\ da 
S F ^(b, k, n) = c F (a, b) T2 h { j h { j , 

with 

A £)F (b) = {aeR + ; \c F {a,b)\ > e u ] , 
and where A\(a,b), ^(a, b) are defined by equation A28\) . 
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We can now state our main result: 

Theorem 4.1 The notations are those of Definition ^. ^ Let F G A £t<Ti d an d v £ (0)1/(2 + 
4/cr)). Provided that e > is sufficiently small, the following assertions hold: 

• |cp(a, 6)| > e v only if there exists some I such that 

(38) (a, b) G Z e = {(a, b), |a|V W (6)| - 1| < A} . 

• For eac/i i = {1, • • • , L} and eac/i pair (a, 6) G 2^ /or which |cp(a, 6)| > e", one /ias /or 
i = 1,2 and for some C > 0: 

|Ai(a,6) - dnpi(b)n et ( b ) \ < Ce v . 

• Moreover for any £ G {1, • • • , L}, there exists C > 0, suc/i £/iaf /or any b G M 2 



lim 

5^0 



1 



27rC^ J§i JB e {e»,n,b) 



S 5 f E {b, k, n)dkdn - A e (b)e w(b)n0 iW 



< Ce u . 



where we denote Be(e u ,n,b) = {k G R ; maxj |/cjn — 9 bi (pi(b)n0 e ( b \ | < e^}, and where 
C^p was introduced in ([6]). 



Proof. See Section E3 



□ 



5 Implementation and numerical experiments 

This section aims to illustrate the efficiency of the monogenic synchrosqueezing transform, 
and its potential applications. We will first describe the computation of the discrete MSST 
from a N x N discrete image, then we will show some examples on artificial and real images. 



5.1 The discrete monogenic synchrosqueezing transform 

The first step consists in the discretization of the variables, therefore one builds discrete 
sets A, JC, O for the scales a, the normalized frequencies (&i,&2), and the orientations 9 
respectively, the frequency resolution in k being denoted by Note that one usually chooses 
a logarithmic scale for a and k, e.g. A = {2^ nv }j = o.-n a -i, so that we have n v coefficients 
per octave. Then, we compute the discrete monogenic wavelet transform cf(cl, 6), and the 
estimations of the instantaneous frequencies Aj.(a, b) and A 2 (a, b) defined in equation (f2"8j) by 
writing for i = 1,2: 

d bt c f (a,b) = f(x)d bi tp a , b (x) dx, 

the other components d bi CTi 1 f(a,b) and d bi c-]z 2 f{a,b) being computed in the same manner. 
The monogenic synchrosqueezing transform consists in a partial reallocation of the mono- 
genic wavelet transform according to space, frequency and orientation parameters. Given a 
threshold 7, and for all a £ A, (k\, k 2 ) G K? ,6 G O, one simply writes 

(39) S F ^{b,ki,k 2 ,n d ) = > . 

n„ ^— ' a 

(\c F (a,b)\ > 7 
\h - Re(Ai(a, b) n e )\ < 
\k 2 -Re(A 2 (a,b)n e )\ < ^ 
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The term 1( ^ 2 ' comes from the measure ^§ in equation (|37p . as we deal with logarithmic 
scales. Note that if we use a logarithmic scale for the frequencies like in [21], the resolution 
Afc depends on the value k. 

This discrete MSST is relatively easy and cheap to compute, but it is not of great interest 
in practice because of the large number of variables which prevents both easy representation 
and fast processing (i.e. decomposition, demodulation). To decrease this number of variables, 
one can remark that if a superposition of IMMFs satisfies the separation condition of equation 
(fM|) . then for any I > £' it also satisfies |V<p^(x)| > |V(/?^(x)|, so that the IMMFs are also 
separated in terms of the norm of the frequency vector. This leads us to define the following 
discrete isotropic MSST for all a G A and k G fC: 

(40) S F>1 (b,k)- l ° s{2) ^ CfM) 



a£A S.t. 




— 2 



V|Ai(a,6)| 2 + |A 2 (a,6)|^ 

Intuitively, Sf^Q), k p ) contains all the coefficients Ci?(a,6) whose "instantaneous isotropic 
frequency" y / \A\(a, b)\ 2 + | A2 (a, 6) | 2 is around the value k. The latter isotropic MSST can 
be particularly useful when visualizing the time-frequency distribution, as it only depends 
on three scalar parameters (bi,b2,k), that is by far easier to represent. Hence, this version 
will be used in the experiments hereafter. Additionally we will use the Morlet wavelet with 
parameters \x > and a > 0, defined in the Fourier domain as follows 

(41) ^ i(T (0 = exp(-7rV(|e|-^ 2 ). 
5.2 Representing a synthetic 3-components signal 

This section illustrates the interest in using the MSST to represent multicomponent signals 
in 2 dimensions. Let us first define the synthetic 3-components test-signal / = f± + f 2 + f'3 
with 

( frfay) = e- 10 « a; - a5 ) 2+ ^- ' 5 )- 2 ))sin(10^(x 2 + y 2 + 2(x + 0.2y)) 

(42) { f 2 (x,y) = 1.2sin(40vr(x + y)) 

[ / 3 ( x , y ) = cos(2vr(70x + 20x 2 + 50y - 20y 2 - 41xy)) 

We compute both its monogenic isotropic wavelet transform cp and its isotropic MSST Sf- 
We then aim at representing |cf| (resp. |»Sf|), which depends on the three scalar parameters 
x, y and a (resp. k). We will here propose two distinct visualizations: either a 3-dimensional 
scatter "density" plot, or a 2-dimensional representation for x or y fixed, which looks similar 
to a ID continuous wavelet representation. The following Figure Q] shows these different 
representations for the wavelet and the MSST of the signal (|42p . We see how the MSST 
sharpens the time-scale representation, making it more readable. 
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5.3 Decomposition and Demodulation 

Once the MSST has been computed, the key point to separate and reconstruct the modes is 
to identify the instantaneous (anisotropic) frequency of each mode \Vipi(b)\ at each point b. 
The t th IMMF Fi> is then approximated by F e (b) as follows: 

(43) F t (b)= S Fn {b,k). 

ke(b)—K<k<k((b)+K 

where ke(b) is an appropriate estimate of |V(/?£(6)|. Obtaining this estimate is an involved 
problem which will not be discussed here, so we simply use the approach used in [21] based on 
a greedy algorithm. Our approximation Ff(b) of the I th IMMF is computed by summing the 
coefficients of the synchrosqueezed transform in the vicinity of this instantaneous frequency, 
to regularize the solution. In practice the width k of the window is set to To evaluate 

the method, we compare the mode F% of the test-signal defined in equation (|42|) to the 
corresponding mode Fj> obtained from the discrete MSST after extraction, by computing the 
normalized Mean-Squared Error (MSE): 

(44) MSE(F e ) = " h ~ f '" 



\\m 

Each reconstructed IMMF is displayed on Figured where one also provides the corresponding 
MSE. It is clear that the modes are reconstructed with very high accuracy, although our 
method for ridge extraction is quite simple. 

Finally, we aim here at showing that the MSST remains efficient on real images or tex- 
tures, even though they can not be considered as multicomponent signals. The following 
example shows the image Lenna where we artificially added an oscillating pattern. After a 
2D synchrosqueezing transform one is able to extract the oscillating component to roughly 
reconstruct the original image. Note that this kind of image was first introduced in |27| to 
illustrate the discrete monogenic wavelet transform. 
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(a) (b) (c) 

Figure 2: The extracted and reconstructed modes 1, 2 and 3 for the test signal of equation 
(|42p . The normalized MSE are respectively 0.03, 0.06 and 0.05. We used the Morlet Wavelet 
ifi^o- defined in (|4"Tj) with a = 2 and \i = 1, and removed the border (| at each side for each 
dimension) to avoid border effects. 




(a) (b) (c) 

Figure 3: Illustration of possible applications on real images, (a): Input image, Lenna plus 
an oscillating pattern, (b): the extracted mode by the 2D MSST. The parameters are the 
same as for Figure [2j The MSE is 0.12. (c): The residual approximates the original Lenna 
image. 

6 Proofs 

6.1 Proof of Proposition 13.11 

We prove Proposition 13.11 into several steps. 



6.1.1 A preliminary result 

We first need a preliminary result: 

Lemma 6.1 Assume that F is an IMMF with accuracy e of the form t!3\) . The following 
estimates hold for any y,h G M 2 : 

(45) \A(y + h)- A(y)\ < e (\h\\V<p{y)\ + \h\ 2 ^ 




(46) \ cos e{y + h)- C o S e{y)\<e[\h\\V^y)\+\h\ 2 ^j , 
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(47) |sin% + /0-sin%)| < e (\h\\V<p(y)\ + H'y) 



(48) \V<p(y + h)- V V {y)\ < e (\h\\V V {y)\ + \h\ 2 ^j 



where M is the constant defined in fTl 
Proof. We only prove inequality ()45|) . the others follows in the same way. Observe that: 

A(y + h)-A(y)= [ V A(y + th) ■ h dt 
Jo 



By assumption (|16p . VA is slowly varying with respect to then: 

\A(y + h)- A(y)\ < e\h\ I \V<p{y + th)\dt 

J 

Applying the usual mean value theorem, we deduce that for any t: 
\V<p(y + th)\ < \V<p(y)\ + \V<p(y + th)-V<p(y)\ 

< \V<p(y)\ + t\h\ max ( sup \d 2 ip{y)\ ) 

< |V^(i/)|+Mt|fe|, 

applying Assumption (fl5j) on second order partial derivatives of (p. To sum up, one has 
\A(y + h)- A(y)\ < e\h\ J \V<p(y)\dt + eM\h\ 2 J tdt < e (\h\\V<p(y)\ + M^-\ , 
which is the required result. □ 

6.1.2 Proof of Proposition I3TT1 

We now prove Proposition 13.11 

Proof of point ([!]) of Proposition 13.11 

By definition of the wavelet coefficients of F, one has 



c F (a,b)= [ A(x)e^ ne Ma- 1 ip(- — dx 
Jm.2 V a J 



We now split this sum into three terms 

c F {a,b) = a- l A(b)j R2 B^ n ^i>(^)dx 

(49) +a~ 1 A{h) J R2 \ e ^ n <>M - ^ x ) n m j ^ [*=±) dx 

+CL- 1 J R2 [A(x) - A(b)] e^W0 (2^) dx . 

Observe now that 

(50) (p(x) = (p(b) + V<p(b) -(x-b)+ [ [V<p(b + t(x - b)) - Vtp(b)] ■ (x - b)dt . 
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We set u = and use equation with k = V(p(b). We then deduce that: 

(51) a' 1 I e Mb)+Vrtb)ix-b))n e(b) ^ fEz]L) dx = e <p(b)n e(b) a ^( aV ^(&)) . 

Jm. 2 V a ) 

Combining equations ([75|) . ([50]) and (f5"Tj) implies: 

c F (a,6)-a^(aV^(6)) e ^ n °m 

= A(b) f e Mb)+V<p(b)-(x-b))n g{b) f e n g(b) fi[Vip(b+t(x-b))-Vip(b)]-(x-b)dt _ A Q -l^ 



x — b 



+A(b) 



^tp(x)n g(x) _ <p(x)ng( b) 



a~ x ip 



x — b 



dx 



+ / [A(x) - A(6)]e v(:c)n ^)a~V 



x — b 



dx . 



Hence 



c F (a,b) - a$(aVip{b)) \A(b) e^ (6 H(&) 

< a _1 yl(6) / e [/ 1 [V V (6+t( a; -fe))-V^(6)]-( a; -fe)dt]n e(i , ) _ 1 

Jr 2 



V' 



x — b 



a 



dx 



+a- 1 / \A(x)-A(b)\ 



x — b 



dx . 



dx 



We now give an upper bound of each of these three terms. We bound the first one using 
mean value theorem and the triangular inequality: 



a^A(b) 



^([V<p(b+t(x-b))-V<p(b)Hx-b))dt]ne (b) 



1 



x — b 



\(y<p(b + t(x - &)) - Vp(&)) • (x - b)\ dt 



We now use inequality (jl8"|) with y = b and h = t(x — b). We then obtain: 



a -1 A(6) 
< ea _1 A(&) 



|(V^(6 + t(x - b)) - V(p(b)) -(x-b)\dt 

)dt 



■4> 



\\t\\x - b\ 2 \V<p(b)\ + M\t\ ^ x n 6|S 



dx 

x — b 
a 

x — b 
a 

x — b 
a 



dx . 

dx 
dx 



x - b\ 2 ,„ „ Jx - b\ 3 



6 



x — b 



dx 



where we set u = s — in the last integral. 

Let us now give an upper bound of the second term. Since 

^(x)n e(x) _ ^(x)n e(b) = sin{lp{x)){ne(x) _ 
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one has: 



jp{x)n e(x) _ ip(x)ng( b ) 



< a 1 A(b) / \n e{x) - n e[b) \ 



x — b 



x — b 



dx 



dx 



1 A{b)V2 [\x-b\\V^{b)\ + \x-b\ 2 —\ 

< ea 2 A(b)V2 [ (\u\\V<p(b)\ + a\u\ 2 — J Mu)\du 

JlR2 V 2 / 

< eA(b)V2 (a 2 \Vy{b)\h + a 3 ^I 2 



x — b 



dx 



(53) 



where we have used inequalities (|46p and (|47p with y = b and h = x — b for the third inequality 
and set u = (x — b)/a for the fourth one. 

Finally, using the same approach we prove that the last term is bounded by, using (I45p : 



(54) 



\A(x)-A(b)\ 



x — b 



M 

dx < e ( a 2 \V(f(b)\h + a 3 —I 2 



Combining inequalities ([52l |53"1 [54^) leads to: 

Ka,6) - a4>(aV<p{b))A(b)e^ n ew | < ea 2 A(6) (Y| |Vp(6)| J 2 + yMJ 3 ) + v 72 f|V^(6)|Ji + 

+ ea 2 ^|V^(6)|Ji+ayJ 2 ^ 



which gives (|26|) . 



x — b 



dx . 



Proof of point ([TT]) of Proposition 13.11 

Since A G L°°(M 2 ) and ip G W 1 ' 1 ^ 2 ), one has: 

a bi c F (a,6) = - / ^(xle^W^Wa- 2 ^^ 

As below we split the sum into three terms 

d bt c F (a,b) = - a - 2 A(b)J R2 e^ n ^d x ^{^)dx 

(55) -a~ 2 A{b) f R2 (e^GOM*) - e*w n «( 6 )) <9 Xi V> (^) dx 

-a~ 2 f R2 [A(x) - A{b)} ert'fr'Wd^ dx 

We again use (f50|) and equation ([24"]) with fc = V(p(b) to write: 

(56) - a" 2 y e Mb)+v V (b)ix-b))n e(b) f^-^\ dx = e ^ ne ^ d bi <p{b)n m a^(aV<^(6)) 
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Hence we deduce that 

d bi c F (^b) - d bi cp(b)n e(b) (a${aV<p(b)j) (A(b)e*® n m 



< a- 2 A(b) 

+ a - 2 A(b)V2 [ 



(fo[Vv(b+t(x-b))~\7<p(b)]-(x-b)dt)n e(b) _ 1 

x — b 



x — b 



dx 



\ n e(x) - n ${b)\ 



dx 



-2 



\A{x)-A(b)\ 



x — b 



dx 



A similar approach than in the proof of point (JTJ) leads to the following upper bound: 

d hi c F (^b) ~ d b .Mb)n e{b) (a^aV^ft))) (A(b)^ n ^) 
< eA{h) (^\ Vlp{b) \^ + ^ M A 

+eA(b)V2 (a|Vv9(6)|/;+a 2 ^/^ 
+e (a\Vtp{b)\l' 1 + a 2 ^-I^ 
which leads to the estimate (Jn]) of Proposition 13.14 with: 

\R 2 (a, 6)| < A(b) C| |V^(6)| l' 2 + ^mA + (flA{b) + l) (\Vip(b)\l[ + a™l' 2 



6.2 Proof of Theorem 

We can now prove Theorem 13.21 By definition for i = 1, 2: 

Ai(a,b) = (d bi c F (a,b)) (cFia^b))' 1 . 

Set B = a^(aVip(b)) (A(b)e v ^ n ^. Proposition O implies that: 

(57) \c F (a,b) -B\ < ea 2 \Ri(a,b)\ , 
and 

(58) \d bi c F (a, b) - d bi <p(b)ng(b)B\ < ea \R 2 (a, b)\ , 

where Ri(a,b) and R 2 (a,b) satisfy inequalities ([26]) and ([2?]) respectively. Then, for i = 1,2: 
ki(a,b)-d bi ip{b)n e[b) = [d bi c F (a, b) - d bi ip{b)n eib) B + d bi ip{b)n e(b) (B - c F (a, b))] [cp^b)}^ 1 
Using ([57|) and ([58]) implies that: 

\ki(a,b) - d bi tp(b)n e{b) \ < e[a\R 2 (a,b)\ + a 2 \Ri(a,b)\\d bi tp{b)\] |cF(a,6)| _1 . 
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By assumption \c F {a, b)\ > e v for all (a, b), which leads to (|29p . The last estimate (|30p follows 
from the fact that for any given ao > 

sup (o|# 2 (a,6)| +o 2 |i?i(a,6)||V^(6)|) < oo . 

(a,b)e(0,oo)xK 

Then, since f G (0, 1/2), there exists some £o > depending on ao such that for any < e < £o 

a\R 2 (a,b)\ + a 2 \Ri(a,b)\\V<p(b)\ < e 2v ~ x , 

which is equivalent to 

e x - v (a|i? 2 (a,6)| +a 2 |i? 1 (a,6)||V V 9(6)|) < e u . 
This last inequality and inequality ([29]) clearly implies inequality (f30|) . 

6.3 Proof of Proposition 13.31 

The proof of Proposition 13.31 relies on Theorem 13.21 and on the following lemma: 

Lemma 6.2 Let F an IMMF with accuracy e > 0. Provided that e > is sufficiently small, 
the sign of Re(db 1 c F (a, b)db 2 c F (a, b)) is this of db 1 ip(b)db 2 (p(b) . 

Proof. Let us first observe that Theorem 13.21 implies that for all (a, b) under consideration: 

Ax (a, b) = d b M b ) n e(b) +0{e u ) 
A 2 (a,b) = d b2 ip{b)n e(b) +0(e v ) 

Then 



(59) Ax (a, 6)A 2 (a, b) = ^(0)^(6) + 0(e" 



Hence, for e > sufficiently small, the sign of Ai(a, 6)A 2 (a, b) is this of db 1 (p(b)db 2 <p(b). 



To get the required conclusion, we now relate Ai(a, 6)A 2 (a, b) and db l CF(a,b)db 2 CF(a,b). 
Let us first remark that the definition of the two Clifford vectors Ai(a, b), A 2 (a, b) from equa- 
tion ([28]) implies that : 



(60) Re ( Ai(a, 6)A 2 (a, b) ) = Re (db l c F (a, b) (c F (a, b)) l db 2 c F (a,b) (c F (a,b)) 1 



We use now the fact that db i c F (a,b) and c F (a,b) are both Clifford vectors. We then apply 
equality ([83]) with q = db 2 c F (a,b) and q' = (c F (a, Hence we get that 



d b2 c F (a,b) (c F (a,b)) 1 = (cf(o, &)) 1 x d b2 c F (a,b) . 
Combining this last equality with (|60p leads to 

(61) Re ( Ai(a, 6)A 2 (a, 6) j = Re (c^c^a, 6) db 2 c F (a, b) ) |cf(o, 6)|~ 2 



Taking into account (|6T|) , equation [59] reads 
Re ( 9^0^(0, 6) db 2 c F (a,b) 
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which implies that Re ydb 1 cp(a, b) db 2 CF{a,b)j and db 1 ip(b)db 2 tp(b) have the same sign for 
e > sufficiently small. □ 

We can now prove Proposition 13.31 
Proof of Proposition 13.31 

Theorem 13.21 directly implies that for any i = 1,2 and (a, b) under consideration: 

||Ai(a,6)|- |^(6)|| <£ v . 

Since 9^(^(6) is assumed to be always positive, we get ([3T]) . Since db 2 (p(b) and db 1 ^p(b)db 2 ^p(b) 
have the same sign, Lemma 16.21 implies (|32p . 



6.4 Proof of Theorem 14.11 

The proof of Theorem 14.11 consists in several steps. 

Lemma 6.3 Let F be a function belonging to A e ,d- For any (a,b), there can be at most one 
£ 6 {1, • • • , L] such that 

(62) |a|V<^0)| - 1| < A . 

Remark 6.1 Condition Ii6ty) is a necessary condition for having ?Jj(aV(pi(b)) ^ 0. Lemma \6.'3\ 
then states there is at most one £ such that ip(aV(pe{b)) ^ 0. Then the two following sums, 
involved in the estimation o/c^(a,6) and d^CF^a, b): 

L 

J^Mb) e w(J) "W a$(aVip e (b)) , 
t=\ 

and 

L 

J2Mb) e^ n °cW a^{aV^{b))d m {b)n ei{b) , 

e=i 

contain at most one term. 



Proof. We follow the same line as in |21| . Assume that there exists some (^1,^2) such 
that ([62]) is satisfied with £\ ^ £2- One may suppose that l\ < £2- One then has for j = 1,2, 



a~ 2 (l-A) 2 < 
which directly implies 



9<Pli 0) 



+ 



dx 2 



< a- 2 (l + A) 2 



+ 



dxo 



&Ph (b) 
dx\ 



+ 



9^2 (b) 
dxo 



> 2a~ 2 (l - Af 



and 



&Ph (b) 
dx\ 



+ 



8x2 



9ft! (b) 
dx\ 



Agfa (ft) 

8X2 



< 4a" 2 A . 
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Further, by assumption (|34|) . for any i = 1,2 one has 



dip e2 (b) 



dxi 



d<Pti 0) 



dxi 



> 



d<pt a (b) 



dxi 



dtp i2 -i(b) 



dxi 



> d 



d(p h (b) 



dxi 



+ 



d<ft 2 -i(b) 



dxi 



where we applied assumption ([35]) : 

1 1 0*0 1 - 1^^2-1(^)11 > \dx^e 2 (x)n ee{x) - d Xi (pi 2 -i(x)ne e ,( x )\ > d [\d Xi <pe 2 (x)\ + \d Xi (pe 2 
The usual inequality (u + v) 2 > u 2 + v 2 , valid for any u, v > implies that 



d(fe 2 (b) 



dxi 



dxi 



> d 



d(p h (b) 



dx; 



dx. 



Summing over i = 1, 2 gets: 

4a -2 A > 2da _2 (l - A) 2 > 2da~ 2 (l - 2A) 

that is 



A > dj [2(1 + <£)] 



which leads to a contradiction. 



□ 



Lemma 6.4 Let F be a superposition of IMMF with accuracy e and separation d of the 
form For all (a, b) £R*x M 2 one has 



c F (a,b) -aJ2 ${aV<pi(b)) (A e {b) e ^ (6)n W 



i'=i 



< ea 2 ri(a,6) 



with 



(63) 



ri(o, b) = h Eii(V2^(6) + 1) |V^(6)| +af [Eti ^(&) (|V W (6)| + Af*) 



(T a freing introduced in 1121]) ). In particular 
• If for some ^ G {l,--- , L] 

(64) \a\V(pe (b)\ - 1| < A . 

f/ien 



(65) 



(a, 6) - a^(aV^ (6))^ (6) e ^ (b)n % (6) < ea^M) 



If for all £ 6 {1, • • • , L}, Condition fl&^[ ) is not satisfied then 
(66) |c F (a,6)| < ea 2 ri(a,6) . 
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Proof. By assumption F is of the form 

L 



e <Pe n i 



=1 



where for each £, Fj> = A(,e tptni is an IMMF of accuracy e. Proposition 13.11 applied successively 
to F\ , ■ ■ ■ ,Fl and the linearity of the continuous wavelet transform imply: 



c F (a,b) - a^2i>(aV<p£(b))Ai(b) e 



<Pt(b)n g (6 ) 



< ea 2 ri(a,6) . 



The rest of the proof follows from Remark 16.11 which states that under condition (|64p , the 
following sum 

L 

Y^Mb) e^ (6 Kw atj>(aV(pt(b)) , 

reduces to A £o (b) e w ° (fc)n V b) a$(aV(p £o (b)), and if condition (JMI is never satisfied. □ 

Lemma 6.5 Let F be a superposition of IMMF with accuracy £ and separation d. For all 
a, b such that 

(67) |a|V^(6)| — 1| < A , 
one has, for any i = 1,2: 

(68) d bt c F (a, b) - ad m (b)n ei{b) ^{aV W (b)) (A e (b)e n{b)n ^) < eaT 2 (a, b) , 



with 
(69) 



T 2 (a, 6) = l'x Eti I V^(6)|( V2^(6) + 1) + a4f ^ti {M t + A e (b)(V2M e + |V^(6)|)) 
(I' a being introduced in IF?])). 



Proof. As in the proof of Lemma 16.4] F is of the form F = Y2i=i A.t^ ptnt > where each 
.F) = Aie Vtni is an IMMF of accuracy e. Proposition 13.11 applied successively to F\, ■ ■ ■ ,Fl 
and the linearity of the continuous wavelet transform imply: 

L 

d bi c F (a,b) - a]T dm'(b)n eAb) ii(aVMb)) (if(4)e"'' (^) > (i, ) < eaT 2 {a,b) . 
i'=i 

Remark 16. II states that under condition (I64p . the following sum 

L 

^ $M*K,(6)#*V^(&)) (^(6)e w(6)n V< 
£'=i 

reduces to d m {b)n dl(b) i>{aV W {b)) (^(&)e<«< 6 Kw) . □ 
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Lemma 6.6 Let F be a superposition of IMMF with accuracy £ > and separation d. Let 
I G {1, • • • , L} . For all a,b such that 

(70) \a\V<p e (b)\-l\<A, 
one has, for any i = 1,2 

(71) \Ha,b) - d bi w{b)n 9e(b) \ < ae 1 -" {r 2 {a,b) + a|V^ (6)|ri(a,6)) . 

Proof. The proof follows the same lines than that of Theorem 13. 2^ replacing the two esti- 
mates (|57p and (|58p on the wavelet coefficients and their partial derivatives, by (|65p and (|68p 
respectively. □ 

We can now prove Theorem 14.11 
Proof of Theorem 14.11 

In the proof of Theorem 14.1} the following notations will be needed: 

1 - A „ _ 1 + A 

^min " 



max te{1 ... ;L} sup b€K 2(|V^(6)|)' max min te{lj ... jL} inf 6eK2 (|V^(6)|) ' 

By assumptions on the phases (pg of the modes of F, a m \ n and a max are both finite and positive. 
For any t G {1, • • • ,L}, e > 0, b G R 2 and n G S 1 , one also defines: 

Bi(e, n, b) = {k G M 2 ; max - d b .(p e (b)ne e (p)\ < e} . 
We now show Theorem 14.11 We set e = e v and fix t G {1, • • • , L}. Let us first prove that 

= 2 ^ fA £ , F (b)n{a-,m aXt \A l (a,b)-d bz Mb)no e(b) \<^ a ~ C H a > & ) da • 

By definition of Sp", we have 

Sp U £ (b, k, n)dkdn 

S 1 JB e (e,n,b) 



JS 1 JB e (e,n,b) \JA £tF (b) fj[ V * / « / 



We first use Fubini Theorem to interchange the order of summations in this integral. Remark 
that: 



Js 1 JBe{e,n,b) J A E<F {b) fj[ \ ° J 

< [ a-\ F {aM I [ r 2 ff h ( *j-MU«,b) n) \ 
JA EF (b) JS 1 JR 2 ~_i V ° J 



^-dfcdn 
a z 



'A s , F (b) 

le,F(6) 



/ a 2 |cf(o, b)\da < + oo , 
J A, F tb) 



25 



where the last inequality is coming from the fact that F is a finite superposition of IMMFs 
Fg, whose wavelet coefficients satisfy (|19|) and (|20p (see Remark 13.11 for more details). Then 
by Fubini Theorem: 



F,e 

S 1 JB e (e,n,b) 



Sp l ' (b, k, n)dkdn 



[ a-*c F {a,b)( [ [ 6-ZT\h( ki - Rei Y a,b)n) )dkdn)da. 



2 

A : ., i.V) V-/:' 1 J6 ; r. „.h) " 

Since the integrand is bounded by: 

a~ 2 c F (a,b) H jT (?wfc) ^ 2 n^( ^" Re( t (a ' b) dMn ) ^ 2vr||/ l ||i 1(R) a- 2 |cHa,6)| 
which belongs to L l (A £jF (b)), the Dominant Convergence Theorem applies. Hence, 
lim / / S^(6, /c, n)d£;dn 

8->0J§i JB e (e,n,b) 

= [ a"V(a,6) (lim / / 6^ f\ h ( ^^^) dkdn) da 

JA e , F (b) \^J § i J Bt (e,n,b) f = i V d J ) 

2 



a- 2 c F (a, b) ( lim / / «T 2 TT J» ( ~ w) ) dfcdn ] da 

A £ , F (6) y^Ojgi JB e (£,n,b)nB'(SR,n,b) fj[ V ° J J 

where we set B'(5R,n,b) = {k £ M. 2 ; maxj |/cjn — Aj(a, 6)| < with i? the bound of the 
support of h. Remark now that, if maxj |Aj(a,6) — db i <pi(b)riQ e ^\ > e, then for 5 sufficiently 
small 

B e (e,n,b)nB'(SR,n,b) = , 

and the limit above is equal to 0. On the other side, if maxj |Aj(a, b) — d^cpiQijng^l < e, for 
5 sufficiently small: 

B e (e, n, b) n B'(SR, n, b) = B'(SR, n, b) , 

and 

lim / / s - 2 fr h (h Z MHa 1 b) n)\ ^ 

J§i JB'(5R,n,b) fJi \ 6 / 

2 



= 2vr, 

where the last equality comes from the assumption L /i = 1. One then deduces equality (f72|) 
Observe now that for any e > sufficiently small there exists some ao > a m ax such that 

(73) 2vr a Q a < e v and sup e x ~ v a (T 2 (a, b) + a|V^(6)|ri(a, b)) < e v . 

00,1 (a,fe)e(0,a )xR 2 
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Indeed, by definition of Ti(a,b),T2(a,b) given respectively in equations ([65]) and ([69]) . there 
exists Cr > not depending on ao such that 

sup a (r 2 (a, b) + a\Vipt{b)\T\(a, b)) < Crag . 

(a,6)e(0,a )x]R 2 



This inequality implies that to prove (J73J) , it is sufficient to find ao such that 
(74) a max < (2vr \\F\\ B -„ e-») 1/a < a < (C^e 2 ^ 1 ) 1 ^ . 



Since v < 1/(2 + 4/<r), one deduces that g^ 2 " -1 )/ 4 /e~ u ' a tends to oo when e — > 0. Hence, for 
e sufficiently small, there exists some ao > a max satisfying (|74p . From now to the end of the 
proof, we fix such ao- 

We now split the integral in the right-hand side of (|72p into two terms 



lim^o / s i I Be (e,n,b) k > n )dkdn 

( 75 ) = 27F X4 e , J ,(6)n{ae(0,ao);max i |A 4 (a,6)-ft iW (6)n tfi( 6 ) |<6} Ci? ( a ' W 

+ 2w JA e , F (b)n{a>a ;ma Xi |A i (a,6)-a i , < ^(6)n^ (b) |<g} CF ( a ' ? ')# " 



We now deal with each of these two terms successively. Recall that since by assumption F is 
a superposition of IMMFs, then inequality (|19p holds for some a > (see Remark 13. ip . We 
then get that for any ao satisfying ([73]) : 



27r X4 e , F (6)n{o>ao;ma Xi [A i (a,6)-a i) . w (6)n 9 . (i , ) |<g} c H a , & )^ 



(76) < 2vr/+ oo | CF (a,6)|^<2 7 r||F|| B - :i a -<e. 

Let us now consider the integral 

I n da 

27T / c F (a,b)— , 

J A e>F (b)n{ae(0,ao); max, |Ai(a,6)-S i , i ^(6)n^( i ,) |<e) a 

and prove that under assumption ([75]) 
(77) 

A £jF (6)n{a G (0,a ), maxlA^a^)-^^^)^^^! < e} = A £ , F {b)n{a G (0,a ); |o|V^(6)|-l| < A} . 

Let us first assume that a G A £t p(b)n{a G (0, ao); |a|V^(6)| — 1| < A}. Since |a| V^(6)| — 
1| < A, Lemma 16, 61 and assumption (j73p imply that 

max | A; (a, 6) - d bi ^i{b)n dl{b) \ < ae 1 ^ (T 2 (a,b) + a\Vip e (b)\Fi(a, b)) < e , 

that is a G {a G (0,a ), max; |Aj(a,6) - ^^(6)n^(6)| < e}> 

Conversely, let us assume that a G A £) i?(?)) H {a G (0, ao); maxj |Aj(a, b) — ^^(6)^(6)1 < 
e}. Since a G j4 £j j7(6), |ci?(a,6)| > e. Then by Lemma E3] there exists some t! G {1, ••• ,L} 
such that |a|V^/(6)| — 1| < A, otherwise assumption ([73]) and equation ([66]) would lead to a 
contradiction. In addition, remark that if £' ^ £, one has 

max|Aj(a,6)-<9 6 ^(&)rt^( 6 )| > \d bi ipi(b)n ei{b) -d bi ip e (b)n 8e , {b) \-max |Aj(a, 6) — (6)n^, (6) | 
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As above, Lemma 16.61 and assumption (|73p imply that 



max | Aj (a, b) - d bi ip v (b)n 6t ,{b)\ < £ ■ 



Hence, we get that 



max|Ai(a,6) - d bi tp e (b)n ei{b) \ > d [\d b .ip £ (b)\ + \d bt (pe(b)\] - e 



> e 



provided that e is sufficiently small, namely that: 



(78) e< [-[\d bi w{b)\ + \d bi w 



i.e. d [|^(6)l + \d bi ip e/ (b)\] > 2e, W,^,Vi 



Hence, necessarily, I = £' and then a G A £ ^(b) n {a G (0, ao), |a|V^(6)| — 1| < A}. 
In addition, one has 



2ir I CF(a,b)a 2 da 

' J 4 e , F (6)n{ae(0,a ),|a|V^(6)|-l|<A} 



2tt 



{ae(0,ao), |a|V<^(6)| 



Ci?(a, 6)a 2 da — 2ir I 
-1|<A} ./{, 



ae(0,a ),|a|V^(b)j-l|<A}\A EiF (6) 



cj?(a, 6)a 



with |/| < e. Using this last equation and equations (|75|) and (jT6j) . we deduce that 



li m J_ f f S S / e {b,k,n)dkdn - yl £ (6)e^ (b)n <W 
5^0 (7^ y§i JB e (e,n,b) 



< 



2tt 



+ 



C^. V{a6(0,o ),|a|V^(6)l-l|<A} 
2vr f 

y{ae(0,ao),|a|V^(6)|-l]<A}\ J 4 e]F (6) 



c F (a, 6)a~ 2 da - A e (b)e Mb)n<> ^ 
cp(a, b)a~ 2 da + e . 



Observe that the domain {a G (0, ao), |a|V(/?£(6)| — 1| < A} is bounded from below by o n 
which only depends on F. In addition, if a G" A s f(6), inequality ([66]) is satisfied. Hence 
(79) 



/■fa 



{ae(0,a ), \a\V(pe{b)\-l\<A}\A s ,F(b) 



CF(a,b)a 2 da 



< e 

< e 

< e 



■CL [ su Pae(o,a )a 2r iM) 
su Pag(o,a ) a 2 ri(a,6) 

^min 

SU Pae(0,a )« 2r i( a ^) 



a 2 da 



l-u 



< Ce, 

where C = (a m ; n inf fceK 2 |V^(6)|) _1 . The last inequality comes from assumption ([73]) valid 
for £ > sufficiently small. This last inequality and inequality (I65p proved in Lemma 16.41 
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then imply that 



< 



lim 

2vr r 

Cij} J{ae(0,a ),\a\V<p e (b)\-l\<A} 

2vr 



+ 



J{a£(0,a ), |a|V</^(ft)|-l|<A} 



S^(b, k, n)dkdn - A e (b)e w{b)n W 

a^(aV^(6))^(6)e^ (6)n <w) a~ 2 da - A t {x)e^ x)ne i^) 
[ea 2 Ti(a, b))a~ 2 da 



+ (C + l)e 



To conclude, observe now that the first term of the right-hand side vanishes. Indeed, the 
wavelet is real isotropic and for any b G M 2 , a — > "ifi(aV(pi(b)) is supported in (0, a max ). Since 
do > amax, we then have: 



do 



/ $(aV<p e (b)) 
J{oe(0,oo),|o|V^(6)|-l|<A} a 



il>(aV<p t (b)) 



da _ J_ 

a 2n 



12 ^ 



iei s 



27T ' 



where CU has been defined in ([6]). The second term of the last inequality can be bounded 
using the same approach than in inequality ()79p . We then get Theorem 14.11 



A Quaternionic calculus 

In this appendix, we give a short introduction to the algebra of quaternions H. More details 
can be found for instance in [32J. 

EI is the real 4-D vector space spanned by {l,i,j,k}, i.e. a quaternion is of the form q = 
Qo + Qi 1 + °2j + 93k, where the algebra product is defined by i 2 = j 2 = k 2 = — 1 and ij = — ji = 
k, jk = -kj = i, ki = -ik = j. 

For a given quaternion 5 = 50 + 5ii + ffej + 53k, one defines: 

• its real part: Re(g) = qo, 

• its vectorial part: Vect(g) = (51,52,93) £ 

If Re(5) = 0, 5 = 5ii + q2] + 93k is called a pure quaternion. 

The conjugate of the quaternion 5 is 5 = 50 — 5ii — 52j — 53k, and its norm (or modulus) is 
defined by: 

\q\ = Vtfl = \j ql + q\ + ql + <& 

For any (5,5') G H 2 , one has: 

55' = q' 5 and \qq'\ = \q\ \q'\ . 

Additionally, the quaternions are a division algebra, which means that any (non-zero) quater- 
nion admits a multiplicative inverse given by: 

(80) q- 1 = . 

\q\ 
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Observe that one has: q 1 = -A™ = (q) . The set of unit quaternion {q £ EI ; |g| = 1} will 
be denoted by S 3 . 

The exponential function is defined on the quaternion algebra as 

+ °° q e 

exp :g ^e g = ^-, 

which converges since exp(|g|) converges. Using the exponential map, each unit quaternion 
of S 3 (\q\ = 1) such that Vect{q) ^ can be written as: 

Vect(q) 

(81) q = (cos ip+n sin ip) = ^ pn with n = — -, cosp = Re(q) and sin ip = \ Vect(q)\ 

\Vect{q)\ 

which is an extension of the complex exponential. Notice that the usual property e^e u = e^ +u 
is no more satisfied in general. The polar form of any quaternion q is then given by: 

(82) q= \q\ (cos (p + n sine/?) = \q\ e^ 71 

where n is a pure unit quaternion: n = ai + 6j + ck with a 2 + b 2 + c 2 = 1, and 

If (<p,n) £Rx§ 3 satisfies ®, one says that <p is a scalar argument of q and n a vectorial 

orientation of q. 

A quaternion of the form q = qo + Qi^ + Q2] with (qo, q±, 52) £R 3 (</4 = 0) is called a Clifford 
vector. Observe that if q = ao+aii+a2j and q' = a^+a'^i+a^i with (ao, 01, 02), a 'i> a 2) ^ ^ 3 
then 

(83) ? = 7xg. 

In this case any vectorial orientation of q is of the form: 

n = ai + b] with a 2 + b 2 = 1 . 

Then there exists 9 £ M such that n = cos i + sin 9 j . This implies that q admits the following 
polar form: 

(84) q=\q\ (cos ip + sin tp cos 9 i + sin ip sin j) = \q\ e^ (cos ei + sine J) . 

If (</?, #) £ M 2 satisfies (|84p . one says that ip is a scalar argument of q and a scalar orientation 
of g. Observe that <p and are also Euler angles of the vector \q\\ £t 3 . 

w 

B Hilbert and Riesz Transform 

In this appendix we briefly recall the Hilbert transform and the analytic signal analysis, and 
then proceed with a natural extension in two dimension, the Riesz transform, used to define 
the 2-D monogenic signal in section [TJ 
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B.l The 1-D Hilbert Transform 

We first define the Hilbert transform of a 1-D real valued function. 

Definition B.l Let f G L 2 (R, R). The Hilbert transform of f , denoted as Hf is defined in 
the Fourier domain as follows: for a.e. £ G R, 

ft?(0 = "i sgn(0/(0 • 

Definition B.2 The complex analytic signal associated to f G L 2 (R, R) is t/ien 

F = f + iHf, 

which reads in the Fourier domain as 

F=(l+sgn(0)7. 

We recall below some properties of the Hilbert transform and the analytic signal. 
Proposition B.l The Hilbert transform is an antisymmetric operator on L 2 (R, R): 

<nhj 2 >=-<h,nf 2 > . 

Proposition B.2 1. f and rif are orthogonal 

2. ||F|| 2 = 2||/|| 2 

3. The spectrum of the analytic signal is one sided. 
B.2 2-D Riesz transform 

We begin with the definition of the 2-D Riesz transform of a real valued image. 

Definition B.3 Let f G L 2 (R 2 ). The Riesz transform of f , denoted as IZf is the vector 
valued function: 




where for any i = 1,2, Ttif is defined in Fourier domain as follows: for a.e. £ G R 2 , 

nTf(0 = -i 47(0 • 

We now present the key properties of 7Z (|33j.[27j). The first ones concern the invariance with 
respect to dilations, translations, and the steerability property (relation with the rotations). 

Proposition B.3 The Riesz transform commutes both with the translation, and the dilation 
operator: for any f G L 2 (R 2 ) ; a > and b G R 2 , one has: 

KD a f = D a lZf and KT b f = T b Kf . 
31 



Proposition B.4 The Riesz transform is steerable, that is, for any f E L 2 (M?) one has 

(85) R e (nf)-r d u(R 9 f) - ^_ shld ni{Ref) + cos en 2 {Ref)) ' 

where rg is the rotation matrix defined by 

Proof. We will prove in the Fourier domain that IZ(Rgf) = rgRg{lZf) , which is equivalent 
to equation (185h . Let us first remark that for a.e. (£l 2 , 



(86) (R e f)^) = f( r ^) • 

Using the definition of the Riesz transform in Fourier domain, one has, for a.e. £ £ R 2 : 



nRd)(0 = r e r e l x f-ij-/( r -^ 



r e (Kf) (r-^) 
r e (Wizf) (0 , 



the last relation coming from (|86p . One then deduces the required result. □ 
The Riesz transform is also a unitary and componentwise antisymmetric operator on L 2 (M 2 ): 

Proposition B.5 For any i £ {1,2}, the i-th component of the Riesz transform IZi is an 
antisymmetric operator, namely for all f,g£ L 2 (M. 2 ): 



(87) < Kif, g > L 2 (R 2)= - < /, TZig 

Since TZ 2 + TZ 2 = —Id, it implies in particular that: 

< 7tf,1lg > L 2( R 2 !R 2)=< TZif,TZig >l2( K 2) + < IZ 2 f, 1Z 2 g >l2 (R 2)=< f,g 
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